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The role of the discontinuity of the exchange-correlation potential of density functional theory is 
studied in the context of electron transport and shown to be intimately related to Coulomb blockade. 
By following the time evolution of an interacting nanojunction attached to biased leads, we find 
that, instead of evolving to a steady state, the system reaches a dynamical state characterized 
by correlation-induced current oscillations. Our results establish a dynamical picture of Coulomb 
blockade manifesting itself as a periodic sequence of charging and discharging of the nanostructure. 
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Coulomb Blockade (CB) [JJ is one of the true hall- 
marks of electron-electron interactions in mesoscopic and 
nanoscale physics. In molecular transport, CB is due to 
an electrostatic barrier induced by the electrons in the 
device which prevents further electrons from tunneling 
in [2] . Progress in the theoretical description and experi- 
mental manipulation of CB is expected to foster advances 
in nanoelectronics and quantum information technology 
[2H5]. On the theoretical side, many aspects of CB can 
be understood in a rather simple way [6]. However, to 
achieve quantitative accuracy in real systems an ab ini- 
tio description is required, something which has not been 
accomplished yet. 

In fact, most current ab initio treatments of transport 
are limited to the steady-state regime and based on the 
Landauer formalism combined with Density Functional 
Theory (DFT) [7 . Within this prescription (L+DFT), 
the agreement with experiment is often poor, particularly 
for devices weakly coupled to leads. The approach has 
also been criticized on fundamental theoretical grounds 
(see e.g. [5]). A main reason for its failure has been iden- 
tified in the shortcomings of typical exchange-correlation 
(XC) functionals of static DFT [9], such as, e.g., the lack 
of a derivative discontinuity in popular Local Density Ap- 
proximations (LDA) and Generalized Gradient Approxi- 
mations [TU] • The importance of the derivative disconti- 
nuity in a static picture of CB for finite systems has been 
investigated in Ref. [IT]. 

As transport is an intrinsically non-equilibrium phe- 
nomenon, the scientific community has progressively 
shifted to time-dependent (TD) approaches in recent 
years (see, e.g., Refs. [321 E]). This shift also recog- 
nizes the fact that significant novel physics occurs in the 
time domain. While different approaches to TD trans- 



port have been developed, TDDFT [T3] offers a compu- 
tationally efficient but still in principle exact method[15). 

In relation to CB, it is then quite natural to address 
the following issues: i) what is the connection between 
a TDDFT description of CB in real time and the one of 
standard steady-state approaches? ii) which features an 
approximate XC functional should have to describe CB? 

This Letter takes a first step in answering these ba- 
sic questions. We present here a TDDFT study of CB 
in the time domain for a correlated single-level quan- 
tum dot (QD) weakly coupled to leads. We use a novel 
XC functional ,16] exhibiting the proper derivative dis- 
continuity and find that the assumption that the system 
evolves towards a steady state is not generally justified. In 
contrast, the discontinuity of the potential leads to self- 
sustained oscillations induced by electron correlations, 
with history-dependent frequencies and amplitudes. Our 
results thus reveal dynamical aspects of CB which are 
not accessible by traditional steady-state approaches. 

The system and its TDDFT time evolution. - The 
Hamiltonian for a single-level QD coupled to two semi- 
infinite, non-interacting ID leads is 

H(t) = H QB + H a +H T + H hias (t). (1) 

a=LM 

In Eq. (§, H a = - £ CT ££ x (^4f i« )t A»,* + h.c), and 
Ht = - E a , 5 (Mink c} a o .co CT + h.c.) respectively describe, 
in standard notation, the leads and their coupling to the 
QD. H a and Ht contain only nearest neighbor hopping 
terms: V (in the left, L, and right, R, lead) and 
(to/from the QD). In the following we take both V and 
Vii„k to be positive. The effect of applying a TD bias is 
contained in H hias (t) = Y^ a ,a Yn^i W a (t)cl aj(T Ci a!<7 . Fi- 
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nally, the QD Hamiltonian is 



Hqu — «cxt / , flQo 



(2) 



with fiQu = cjg.coo- the density for electrons with spin a. 
Two parameters describe the QD: the charging energy U 
and a static gate voltage v cx t- For time t < 0, the system 
is in equilibrium; at t > 0, a bias W a (t) is applied. 

Within TDDFT, the many-body Hamiltonian Q is 
mapped onto an effective one-particle Kohn-Sham (KS) 
Hamiltonian. The KS dynamics produces the exact 
TD density, provided we know the exact KS potential 
whose density dependence is non-local in space and time. 
In practice, one has to resort to approximations. As 
shown below, we can already gain significant insight via 
an approximate, adiabatic KS potential. The approx- 
imation we use is based on a ground state XC func- 
tional obtained from an LDA of the non-uniform ID 
Hubbard model, via the Bethe Ansatz (BALDA) [Trj] , 
In the context of TDDFT, an adiabatic version of this 
functional (ABALDA) has already been used to inves- 
tigate the dynamics of Hubbard clusters [T7]- Within 
the ABALDA, the KS Hamiltonian reads 2?ks(*) = 

#QD,KS(*) + Ea=L R H a + H T + with 



(3) 



The KS potential is a functional of the instantaneous 
density no(t) = X}a- n 0crW on the QD. Explicitly, 

wks ["oft)] = v cxt + -Un (t) +v KC [n (t)}. (4) 

The original BALDA functional was devised for a uni- 
form chain with Hubbard interaction U and nearest 
neighbor hopping V: in this case, the ratio U/V is the rel- 
evant parameter in v xc [16 . However, in the CB regime, 
which we wish to address here, the electrons need to be 
largely localized in the device. For our model system, 
this corresponds to i) taking Vu n k significantly smaller 
than V and ii) using U/Vn nk as relevant parameter in 
v xc . Accordingly, we propose the functional form 



v xc [n] = 6(1- n)vg> [n] - 6(n - !)«<<> [2 - n] 



v { < ] [n] = --Un -2V link 



cos 



/7rn\ 
V~2~J 



cos 



(5) 
(6) 



where the parameter /3 is determined by the condition 



2/3 



sin(7r//3) = 4 / dx 



J (x)Ji(x) 



x[l+eMUx/(2V link ))] 



(7) 



with Ji=o t i(x) Bessel functions. For Viink — > 0, w xc [".] — > 
-9(1 - n)Un/2 + 0(n - l)U(2 - n)/2, which correctly 
reproduces the exact XC potential of the isolated QD. 
This v xc also leads to densities (not shown) in good 
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FIG. 1: (Color online) Time evolution of the density, current 
and KS potential for three different biases. In all panels, solid, 
chain and dashed line refer to Wl — 1.3, 1.6, 1.9, respectively. 
Top panel: densities at the QD. The inset shows the density 
at the end of the propagation period. Middle three panels: 
current through the QD (thick lines) and KS potential (thin 
lines). Bottom panel: current five sites away from the QD. 



agreement with the Quantum Monte Carlo data [TS] for 
the Anderson model. At half- filling, d xc [ji] is discontin- 
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4Viink cos(|). 

exact XC potential is certainly discontinuous for the 1- 
D Hubbard model or for the isolated QD. However, for 
the QD connected to non-interacting leads it is reason- 
able to expect a small smoothening of the discontinuity 
(this is supported by preliminary calculations, not shown 
here, of the exact v xc for an Anderson impurity in small 
clusters). We therefore modified the discontinuous (at 
n = 1) t>xcM of Eq. ([5]), with a softened, continuous ver- 
sion v xc (n) = f(n)v£F[n] - (1 - f(n))v ( X c ) [2 - n]. Here, 
f(n) = cxp (( n _\)/ a ) +1 , with a smoothening parameter a. 
The smoothening of the KS potential also has the advan- 
tage of alleviating the numerical difficulties caused by the 
sudden changes of v xc during time propagation. 

To propagate the lead-dot-lead KS system in time, we 
use a recently developed TD algorithm j^Hj for an open 
system of (effectively) noninteracting electrons, as re- 
quired in TDKS. In all simulations below, energies are 
measured in units of V, times in units of V~ x and cur- 
rents in units of |e| V where e is the charge of the carriers. 
The sharp slope near n — 1 in v xc has a profound impact 
on the time evolution of the density on the QD as well 
as on the current through it. 

TD transport in the CB regime.- In Fig. [T] we show 
the time evolution of the density, the current and the 
KS potential with a smoothening parameter a = 10~ 4 
for three bias values and Viink = 0.3, v cxt — 2, U = 2, 
and Fermi energy Ep — 1.5 (as shown below, this choice 
of the parameters corresponds to the CB regime). The 
system is in its ground state for t < 0; the equilibrium 
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density is calculated self-consistently with the vks of Eq. 
Q but in static DFT. At t = 0, a dc bias W L is suddenly 
switched on in the left lead. Remarkably, for this set of 
parameters and within a certain range of biases, the sys- 
tem does not evolve towards a steady state; instead, after 
a transient period, we see self-sustained density oscilla- 
tions around unity (Fig. [Ij top panel) , with an amplitude 
of the order of 10 -3 (inset in the same panel). For the 
lower bias, Wl = 1.3, the density remains below unity 
for most of the period, for Wl =1.6 the density roughly 
oscillates around 1, and for Wl = 1-9 the density ex- 
ceeds unity for most of the period. The oscillations are 
also present in the current (middle three panels), but this 
time with a much larger amplitude relative to its aver- 
age value. The saw-tooth structure is consistent with the 
continuity equation (the density is approximately piece- 
wise parabolic). Note that away from the QD, the os- 
cillations in the current (and in the density) disappear 
(bottom panel). Oscillations are clearly visible also in 
vks (middle three panels). Due to the fairly large jump 
in u xc and the small amplitude of the density oscilla- 
tions, the KS potential is a train of almost rectangular 
pulses. For larger biases the density remains above unity 
for longer times, and the width of the pulses in i>ks ex ~ 
tends in time. The oscillations discussed here are a direct 
consequence of the discontinuity in i> xc . For similar calcu- 
lations (not shown) in the Hartree approximation, where 
Uks is continuous in the density, the system does evolve 
towards a steady state. 

History dependence, oscillations and smoothening.- We 
have calculated the time-evolution of the system as ob- 
tained by switching the bias as WL(t) = Wl sin 2 ( 2T 7r * t - ) 
for t < T sw itch and IMl(£) = Wl — const otherwise. 
The Fourier analysis of no(t) reveals peaks whose posi- 
tion and height depend on T sw i tc h, i-e., on the history of 
the applied bias. We found that for any small but finite 
a, the amplitude of the oscillations approaches zero as 
Tswitch — > oo. We therefore conclude that in this case the 
L+DFT approach gives the same solution as TDDFT for 
sufficiently slow switching [21]. For a = 0, instead, the 
system never reaches a steady state. 

CB regime from the steady-state approach.- The oscil- 
lations just described are distinct from those occurring 
in the presence of single particle bound states (2H |5S] : 
They are induced by electron correlations, i.e., are ab- 
sent in the mean field (Hartree) approximation. An es- 
pecially revealing feature is that at long times the system 
is in a dynamical state of oscillating density and current, 
whose time-average is fairly constant for a large range 
of biases. For further insight, we now use the L+DFT 
approach which, as said above, yields the same solution 
as TDDFT for a + 1 and for adiabatic switching. Using 
non-equilibrium Green's function techniques, the value 
n°° of the steady-state density at the QD is given by the 




FIG. 2: (Color online) Left panel: Graphical solution of 
Eq. (JsJ) for few values of the applied bias Wl- Right panel: 
Steady-state density as function of Wl for a smoothened vks 
(a — 10~ 4 ) and few values of the dot-lead hopping Vnnk- 

self-consistency condition [26] 

/£F + W a J 
— v{u-w a )\G{^ . (8) 

Here = (w - v KS [n°°] - J2 a s (^ ~ W^y 1 is the 

retarded Green function at the QD site, which depends 
on n°° through vks, ^( w ) is the embedding self energy 
of the leads, and the width function T(ui) = — 2Im[E(u;)]. 

We first consider a discontinuous v xc (i.e. a = 0). In 
Fig. [2] left panel, we display the l.h.s. and the r.h.s. of 
Eq. Q as function of n°° for the parameters V = 1, 
Mink = 0.3, U = 2, w cxt = 2, ep = 1.5 and for different 
values of the dc bias Wl ■ These are the parameters used 
in our TD simulations of Fig. [I] Interestingly, there is a 
range of applied bias voltages for which Eq. Q does not 
have a solution. Since Eq. ([8]) is a condition for the QD 
density at the steady-state, it follows that, as a direct 
consequence of the discontinuity in v xc , for some values 
of W a , no steady state exists. On the other hand, for a = 
10~ 4 , as in our TD simulations, Eq. (|8j admits steady 
state solutions; however, from Fig. [I] we know that, in 
general, the system time-evolves towards an oscillating 
regime and that the solution of the L+DFT scheme is a 
solution in TDDFT for adiabatic switching. 

To show that in fact we are in the CB regime, in Fig. [2] 
right panel, we display the steady state density n°° , as 
a function of Wl- As before, u ext = 2, e F = 1.5, U = 2 
and a = 10~ 4 . We see a plateau in the density at n = 1 
developing for a range of Wl ■ The plateau increases with 
decreasing Mink, and in the limit of small Mink its length 
becomes equal to the charging energy U. This is the 
usual CB picture: if the site is occupied by one electron, 
a second electron can only jump in if its energy exceeds 
the charging energy of the QD. For Mink = 0.3, the Wl's 
of Fig[l] correspond to the beginning, middle, and end of 
the density plateau. Our results show that there exists a 
range of biases, starting at a critical value W c , for which 
the decay time of the oscillations becomes infinite. 

Some considerations on v xc and the ABALDA.- Drop- 
ping only the assumption of locality in space, v xc will 
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depend on the instantaneous densities at few, say M, 
sites around the QD. This gives a steady state condition 
similar to Eq. ([8]), in the form of a set of equations for the 
densities at the M sites. For a discontinuous w xc at one 
(or more) of the M sites, these equations may not have a 
solution, similarly to what we found in Fig. [2j The role 
of non-adiabaticity is much harder to assess since it re- 
quires the development of XC functionals with memory, 
e.g., via non-equilibrium Many-Body Perturbation The- 
ory [T31 [STJ HH] or fluid dynamical considerations [Ml [30] . 

We can now provide an answer to the two questions 
posed at the beginning of this Letter. A discontinuous 
or rapidly varying v xc is central to the description of CB 
within L+DFT and TDDFT. For KS potentials with a 
true discontinuity, the steady-state self-consistency con- 
dition in L+DFT cannot always be satisfied. Consider- 
ing a smoothened discontinuity (which, as said earlier, is 
physically more realistic) the L+DFT approach yields a 
clear-cut CB scenario. However, the very same XC po- 
tential used in a TDDFT framework leads to a dynamical 
state with history-dependent, self-sustained density and 
current oscillations: where the static approach gives CB, 
the TD approach gives oscillations. The average of these 
oscillations corresponds to the L+DFT solution. The lat- 
ter is a TDDFT solution only for an adiabatic switch-on. 

To conclude, our results suggest that in a transport 
setup CB manifests itself as an intrinsically TD phe- 
nomenon, i.e., it corresponds to an oscillatory current 
representing the intuitive picture of CB as a sequence of 
charging and discharging of the weakly coupled molecule 
or QD. While our calculations are performed for a model 
consisting of a single correlated QD coupled to non- 
interacting leads, we expect the results to be of gen- 
eral nature: In the continuum-real-space TDDFT de- 
scription, whenever the particle number on the molecule 
or QD crosses an integer, the discontinuity of the time- 
dependent XC potential [3TJ will trigger persisting charge 
and current oscillations which cannot be captured in any 
steady-state approach. 
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